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We investigate the superfluid behavior of a Bose-Einstein condensate of ®Li molecules. In the 
experiment by Weimer et at, Phys. Rev. Lett. 114, 095301 (2015) a condensate is stirred by a 
weak, red-detuned laser beam along a circular path around the trap center. The rate of induced 
heating increases steeply above a velocity Vc, which we define as the critical velocity. Below this 
velocity, the moving beam creates almost no heating. In this paper, we demonstrate a quantitative 
understanding of the critical velocity. Using both numerical and analytical methods, we identify 
the non-zero temperature, the circular motion of the stirrer, and the density profile of the cloud as 
key factors influencing the magnitude of Uc. A direct comparison to the experimental data shows 
excellent agreement. 

PACS numbers: 67.85.-d, 03.75.Kk 


I. INTRODUCTION 

The intriguing phenomenon of superfluidity arises from ^ 
an interplay of quantum statistics and interactions at low ^ 

temperatures. Its defining feature is the stability against ^ 

external perturbations and, closely related, dissipation- 
less flow near an obstacle. However, the superfluid be- ^ 
havior will only be sustained within a certain parameter 2 
regime. If the system is perturbed by a sufficiently strong g 
or sufficiently fast moving perturbation, its dissipation¬ 
less nature will break down. The onset of dissipation 
typically occurs around a velocity which is referred to as 
the critical velocity. 

In a seminal study, Landau estimated the critical veloc¬ 
ity by considering the onset of dissipation due to the ex¬ 
citation of single elementary quasi-particles [T] of the sys¬ 
tem. In order to satisfy both energy and momentum con¬ 
servation, an elementary excitation of energy e(k) with 
momentum fik can only be created above the velocity 
Vc = mini(;[e(k)/(IiA:)], which is the Landau criterion, h 
is the Planck constant and k is the wave vector, with 
k = |k|. When this expression is applied to the Bogoli- 
ubov approximation of the spectrum of a weakly inter¬ 
acting Bose gas, it predicts the sound velocity up as the 
critical velocity. This is the natural scale for the critical 
velocity, if the decay mechanism is due to phononic exci¬ 
tations. In addition to phonons, however, other types of 
excitations can control the critical velocity. In Ref. [5], 
Feynman considered the superfluid flow into a long chan¬ 
nel with a diameter D ^ ^ [3], and estimated a critical 
velocity of Vc ~ h/(rnD) h\{D/^). ^ is the healing length 
and m is the atomic mass. Thus, vortex excitations can 
result in a lower critical velocity than phonon excitations, 
if this velocity scale is smaller than the phonon velocity. 

Several experiments have probed the superfluidity of 
dilute Bose gases via a local perturbation, in particular 
via laser stirring. In Ref. [4] a three dimensional (3D) 
condensate and in Ref. a two dimensional condensate 
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FIG. 1. (Color online) Simulated heating rates for a homoge¬ 
neous condensate at T = 0.07Tc (blue squares), stirred along 
a linear path by a weak stirrer of strength Vo = —0.03 yi. A 
critical velocity of 4.2 mm/s is determined using the fitting 
function given by Eq. The fitted curve is shown by the 
solid line. The Bogoliubov estimate of the sound velocity ub 
is 4.4 mm/s. The insets are snap shots of the density: In 
panel (a) the moving perturbation causes no drag for v < Vc, 
whereas panel (b) demonstrates the formation of sound waves 
leading to heating for v > Vc- The stirrer location is marked 
by a circle of radius a. For the thermal gas with T > Tc 
(red circles) heating increases linearly in v. Here, we use 
Vo = -0.35 At. 


were stirred with a blue-detuned laser beam, moving on 
a linear and circular trajectory, respectively. The break¬ 
down of superfluid flow due to a constriction was probed 
in Ref. [B]: a toroidal Bose-Einstein condensate (BEC) 
was put into quantized rotation, and the decay due to a 
potential barrier, created via a blue-detuned laser, was 
observed. It was demonstrated in Ref. [7], that the on¬ 
set of heating is governed by thermally activated phase 
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slips, due to vortices traversing the barrier region. Re¬ 
cently, superfluidity was observed in an elongated ®Li gas 
oscillating with respect to a ^Li BEC [S] . Here, the onset 
of heating is predicted to occur for a relative velocity that 
equals the sum of the individual sound velocities [3]. In 
addition to these studies, superfluidity in BECs has been 
studied in a number of theoretical works based on the 
Gross-Pitaevskii equation iinHia- Further studies were 
reported in Refs. M- 


As a new, and previously unexplored feature of stir¬ 
ring experiments. Ref. m reported on perturbing a 
cloud of ®Li atoms and molecules, across the BEC-BCS 
cross-over, with a red-detuned laser, resulting in an at¬ 
tractive stirring potential. In this paper, we demonstrate 
a quantitative understanding of the superfluid response 
of the condensate of ®Li molecules. We first demon¬ 
strate that the seemingly minor detail of stirring with 
a red-detuned laser, rather than a blue-detuned one, re¬ 
sults in a qualitatively different scenario for the break¬ 
down of superfluidity. We show that for a blue-detuned 
laser, and for intermediate laser intensity, the primary 
dissipation mechanism consists of shedding vortex-anti¬ 
vortex pairs, rather than fully deconfined vortices, as in 
the Feynman scenario. Interestingly, for a red-detuned 
laser, neither vortices nor vortex-anti-vortex pairs are 
generated. Therefore, the choice of a red-detuned laser, 
with weak or intermediate intensity, compared to the 
mean-field energy, allows to study the phononic dissipa¬ 
tion mechanism with clarity. We identify the condensate 
temperature, the circular motion of the stirrer which was 
used in Ref. m, and the inhomogeneous density of the 
cloud as key parameters that influence the magnitude of 
Vc- These factors reduce the critical velocity below the 
phonon velocity, even though the dissipation mechanism 
is phononic. A comparison of our numerical and analyt¬ 
ical results, and the experimental data of Ref. |15j in 
the BEC regime, gives good agreement, and develops a 
consistent physical picture. 


This paper is organized as follows: In Sec. |^we de¬ 
scribe the simulation method that we use. In Sec. IIIII 
we develop an analytical expression for the heating rate 
within the Bogoliubov approximation, for linear stirring, 
and compare this result to the simulation in the corre¬ 
sponding regime. In Sec. |IV| we study the critical ve¬ 
locity for a repulsive stirring potential and contrast the 
resulting dissipation mechanism to that of an attractive 
stirring potential. In Sec. we develop an analytical 
expression for stirring along a circular path, and com¬ 
pare this result to the simulation. In Sec. |V^we discuss 
the reduction of Vc for finite temperatures. In Sec. m 
we identify the influence of the inhomogeneous density 
distribution of the molecular cloud on Vc- In Sec. |VIII| 
we compare the numerical and analytical results to the 
experiment, and in Sec. IX we conclude. 


II. SIMULATION METHOD 


We simulate the stirring dynamics of the system within 
a c-field formalism. The Hamiltonian of the unperturbed 
system is 


Ho= dr 


+ ^'^Hr)ip\r)^{r)i){r) 
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where "0 ('0^) is the field annihilation (creation) operator, 
and g = Airash^/m is the interaction constant; as is the 
s-wave scattering length. The external potential V (r) de¬ 
scribes the harmonic trap, Utrap(r) = +u!^z^)/2. 

ojr and oJz are the trap frequencies in the radial and trans¬ 
verse directions, respectively, and r = (x^ is the 

radial coordinate. We add the following term to describe 
laser stirring 


Tis (t) = y dr l/(r, t)n(r), (2) 

where V (r, t) is the time-dependent stirring potential, 
h(r) is the density operator at the location r = (x, y, z). 
The stirring potential is a Gaussian along the x and y 
direction with a width cr, and independent of the z direc¬ 
tion: 


with a strength Vq. It is centered around Xs(t) and 
which move either along a linear or a circular path as a 
function of time t. 

To simulate this system numerically, we map it on a 
lattice system, which also introduces a short range cut¬ 
off, see Appendix]^ We represent both the equations of 
motion and the initial state within a c-number represen¬ 
tation, which corresponds to formally replacing the oper¬ 
ators 0 by complex numbers. Furthermore, we approxi¬ 
mate the initial ensemble by a classical ensemble, rather 
than a non-classical phase-space distribution, within a 
grand canonical ensemble of temperature T and chemi¬ 
cal potential g. This can also be seen as the classical limit 
of the truncated Wigner approximation, as discussed and 
applied in Ref. [7], which also corresponds to molecular 
dynamics in the classical wave limit. The initial states 
are generated via a classical Metropolis algorithm. 

Motivated by the experimental setup of Ref. [15], 
we consider a pancake-shaped, three-dimensional con¬ 
densate. The density distribution resembles an oblate 
spheroid. We employ a lattice with 140 x 140 x 11 sites. 
For comparison, we consider a system with homogeneous 
density, without a trapping potential V (r), for which we 
use a lattice with 60 x 60 x 3 sites. For the discretization 
length I we choose 1 /rm. As discussed in Ref. m , I is 
chosen such that the healing length ^ = fi/y/2mgn and 
the thermal de Broglie wavelength, A = ^J2^Th?■ /mkrT 










3 




FIG. 2. (Color online) We show the nnmerically determined heating rate (red squares connected with a dashed line), that is 
induced by stirring with a weak stirring potential of strength Vo ~ —0.03 fi, in comparison to the analytical prediction (green, 
continuous lines). In panel (a), we show the case of a linear stirring motion. In panel (b) we show the heating rate induced by 
a circular stirring motion as a function of the stirring velocity, for a stirring radius of Ro = 10 /rm. In addition to the analytical 
prediction for the continuum case, we also show the heating for a lattice system (green, dotted lines). The inset of panel (b) 
shows the predictions for three different radii, in particular 10 fim, 20 fim, and 30 fim. As the radius is increased, the heating 
rate due to a circular motion slowly recovers the heating rate for a linear motion. Overall, the analytical results capture the 
numerical results well. 


are comparable to, or larger than I, n being the density, 
and kB the Boltzmann constant. In Ref. [15] the con¬ 
densate is stirred using a narrowly focused red-detuned 
laser beam that forms an attractive Gaussian potential, 
whose width is of the order of the healing length For 
simulations of this system we adjust the parameters ac¬ 
cording to the experimental choices, in particular we ad¬ 
just the central column density, the scattering length as, 
the temperature, the stirring time, the stirrer strength 
Vo, and the stirring pattern. Typical trap frequencies 
are ujr ^ 2tt x 30 Hz and Wz « 27r x 550 Hz. For simu¬ 
lations of the homogeneous system we use cr = 1.1/rm, 
n = 0.48 /im“^, and Og = 2180ao, where Oq is the Bohr 
radius. These parameters are in the typical range of the 
experimental parameters. The stirring sequence is the 
following: We linearly turn on the stirring potential over 
10 ms, let it stir the system for 200 ms, and then turn it 
off over 5 ms. This is again inspired by the experimen¬ 
tal choices. We repeat this for various stirring velocities 
V and calculate a thermal ensemble of the total energy 
E = {Hq) using the unperturbed Hamiltonian in Eq. H 
During the stirring process the energy increases linearly. 
From this linear increase, we determine the heating rate 
dE/dt for each v, and for the desired parameter regimes. 
We elaborate on the numerical method of determining 
the heating rate in Appendix [Xj 


III. HOMOGENEOUS BEG WITH LINEAR 
STIRRING 

While the experimental setup displays an interplay of 
features which complicate the interpretation, such as the 
temperature, the stirring pattern and the density inho¬ 
mogeneity, we isolate these features in our analysis to 
identify how they affect the heating rate. We start out 
with the case of a homogeneous condensate at a low tem¬ 
perature, stirred along a linear path with a weak stir¬ 
rer. In Fig. [ijwe show the simulated heating rate as a 
function of the stirring velocity v, for a condensed and 
a thermal system. For the condensed system, the mov¬ 
ing stirrer creates almost no dissipation at low veloci¬ 
ties V. For velocities above the sound velocity ub, where 
= \/9^1 va is the Bogoliubov estimate, sound-wave 
excitations are created and the rate of induced heating 
increases steeply. For these conditions, in particular for 
a temperature much smaller than the mean-field energy 
and for a weak, linearly moving stirrer, we recover the 
estimate Uc ~ t'B- For comparison, we show the heating 
rate of a stirred, thermal gas at a temperature T > Tc in 
Fig. [3 For this regime, we recover a linear dependence, 
corresponding to friction that scales linear in v. 

Throughout this paper, to quantify the critical veloc¬ 
ity, we use the following fitting function of the heating 
rate 


V dt J Fit V 


(4) 


which was discussed in Ref. ini, for V > Vc, with the free 
parameters A, B, and Vc- This fitting procedure gives a 
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in Eq. with a linear stirring trajectory rs{t) = vt, with 
rs(t) = (xs{t),ys(t)) and v = (vx,Vy). In momentnm 
representation, the stirring term is 

=^Vk(iK. (7) 

k 

Vk is the Fourier transform of the stirring potential with 
V = 0, in particular 14 = jA x exp(—A:^(T^/2) 

with A being the system area. 14 (t) is the moving stir¬ 
ring potential, which is I4(t) = 14 exp(ikvt). nj. is the 
Fourier transform of the density operator. We use the 
Bogoliubov approximation to evaluate the heating rate. 
So the Hamiltonian is approximately 
1 2 3 4 i;b 5 .. 

„[mm/s] 'Ho = 2_^fuJkblbk, (8) 
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FIG. 3. (Color online) Simulated heating rates for attractive 
and repulsive stirrers are obtained for various stirrer strengths 
Vo- The negative Vo stands for attractive stirrer, whereas the 
positive Vo for repulsive stirrer. The inset shows the heating 
rates for the repulsive stirrer. 


robust estimate of Vc for both numerical and analytical 
results, even though the analytical form that we find dif¬ 
fers from this expression. For the heating rate of the con¬ 
densed system in Fig. we determine a critical velocity 
of 4.2mm/s, which is comparable to fe = 4.4mm/s. 

In addition to the numerical results in this paper, we 
develop an analytical expression for the heating rate that 
quantitatively describes the simulation results. We derive 
the heating rate using the stirring potential as a pertur¬ 
bation Hamiltonian that we linearize within the Bogoli¬ 
ubov approximation. We extend this to finite tempera¬ 
tures in Sec. |vg where we include the thermal broaden¬ 
ing of the phonon modes. 

In a general setting, not necessarily for a condensate 
or even bosons, we determine the heating rate for a weak 
stirrer perturbatively in the stirring term: 


d{Ho{t)) 

dt 


dh{[Id,j{h),[Hs,iit),Ho]])-.. 


(5) 


where we gave the perturbative expansion to second or¬ 
der. Hq is the Hamiltonian of the unperturbed system, 
such as Eq. Ill and therefore a conserved quantity. /(t) 
is the perturbation in the interaction picture 

= exp{iHot)Hsit) exp{-iHot). (6) 


Hg (t) is the stirring term in the Schrddinger picture, such 
as Eq. The expectation value (...) is taken with re¬ 
gards to a density operator p, i.e. (...) = Tr(p...). If p 
commutes with Hq, as it does for a canonical ensemble, 
the first order term in Eq. [^vanishes. 

We now evaluate this expression for a linear stirring 
motion, also see [TH]. We consider a stirring potential as 


where 4: are the Bogoliubov operators, and huJk is 

the Bogoliubov dispersion, fnok = \/efc(efc -I- 2gnc), where 
Efc is the dispersion of the non-interacting system, and 
ric is the condensate density. The density operator is 
h-k « y/Noiuk + Vk)(p_]^ + &k) with No being the number 
of condensed atoms, u^, Vk are the Bogoliubov functions, 
with {uk + Vk)^ = Ck/ifiuJk)- With this, Hsit) is 

'Hsit) ^ ^ + Vk)ib^^^ + Sk). (9) 

k 


'Hsjit) has the same form with 6k —t 6k exp(—and 
b’'_y. —?> b''_y^exp{iojkt)- With this, we obtain 

dh{[Hsih),[Hsit),Ho]]) (10) 

= + (11) 
^ Wfc — KV 

which at long times approaches 

^ +Vk)‘^No\Vi^\^S{ujk - vk). (12) 

k 

We note that this result can also be obtained by solving 
the dynamical evolution that is created by Eqs. [^and[^ 
exactly, which is 

6k(t) = e-*‘"''‘6k + 7lk(t), (13) 

with 

—21, , /— sinf(a;fe — vk)t/2l 

= —in, + (,.1-vk) 



Therefore the energy at time t compared to the initial 
energy is 


{Hoit)) - {'Horn 


p y] hiOk{uk + Vk)‘^No\V\^\ 

k 

sin^ [(wfc — vk)t/2] 

(Wfc - vk)2 


2 


(15) 
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FIG. 4. (Color online) Analytical heating rate Rn for linear 
stirring, Eq. [^and circular stirring with radius Ro = 10 pm, 
Eq. |33[ at zero temperature, within the Bogoliubov estimate. 
The stirring velocity v is shown in units of ub. 


From this, we recover the heating rate given in Eq. |12| 
It is a sum of Fermi’s Golden rule terms multiplied by 
the Bogoliubov excitation energy. We use {uk + Vk)'^ = 
/{2mLOk) and the explicit form of |Ekp, which results 
in 


R 


SC 


mA 


d^kk^ exp{—k^a^)6{(jJk — vkx), (16) 


where = k^ + ky. We defined the rescaled heating 
rate i?sc = idE/dt)/{NoVQ). Without loss of gener¬ 
ality, we have set the direction of the stirring velocity 
along the x-axis, v = ve^. In addition to this expres¬ 
sion for the continuum system, we also evaluate and 
rescale the expression for the lattice system. We use 
Cfc = 2J— cos(fcjZ)) and numerically evaluate Eq. 
[^for the first Brillouin zone, d is the spatial dimension, 
J = h^/{2mP) is the energy, and I is the discretization 
length. We show the heating rates and their comparison 
to the simulated rates as a function of v in Fig. [^a) . 

In the continuum limit, the heating rate in Eq. |l2| can 
be evaluated explicitly as follows: We write the integral 
in polar coordinates, and integrate out the polar angle. 
This gives 


27rtT^ exp(—/c^cr^) 

= -w / dk — , 

Jo yjv'^k? - ujI 


(17) 


for u > ub, and i?sc = 0 else. The wavevector ku 
corresponds to the maximal phonon momentum that 
can be created by the stirring process, it is given by 
ku = ^“^y/2(u2/u| — 1). We write the heating rate in di¬ 
mensionless form, by defining = h{dE/dt)/{NcyiV^). 
Ncy\ = /A is the number of atoms in a cylinder of 

length Lz and radius cr. With k = k/a, the heating rate 
is 


Rn = 's/S 




dk 


exp(—fc^) 


X-k^l2 


(18) 


with X = (u^/ug — l)(cr^/^^). The integration can be 
performed analytically, giving 

i?„ = 2^Xexp(-X)(/o(X)-/i(X)), (19) 

where Io{x) and Ii {x) are the modified Bessel functions of 
the first kind of order 0 and 1, respectively. This function 
is shown in Eig. The onset of dissipation occurs at 
u = ub, with a steep, linear slope: 

Rn « 2TT{a^ie){vVvl - 1) (20) 

for u/ub > 1. Eor large velocities, v/vb ^ 1, the heating 
rate falls off as 

Rn « v^(^/a)^^^=. (21) 

The heating rate undergoes a maximum at 

<ax«^l(l + 0.79^V^"), (22) 

where it assumes the value (i?n)max ~ 1.65. 

IV. ATTRACTIVE VERSUS REPULSIVE 
STIRRERS 

In the experiment reported in Ref. the stirring 

potential was attractive, in contrast to previous mea¬ 
surements such as Refs. [H [S]. This superficially mi¬ 
nor difference results in a qualitatively different dissipa¬ 
tion mechanism as the stirring potential is increased from 
small magnitudes to the larger ones of the order of the 
mean-field energy. In Eig. we show the heating rate 
as a function of the stirring velocity for increasing stir¬ 
ring potentials, at a low temperature of T = 0.07 T^, and 
for a linear stirring motion. For small stirring potentials, 
compared to the mean-field energy /i = gn, the heating 
rate for attractive and repulsive potentials agree. For 
this limit the response is quadratic in the potential, as 
demonstrated in the previous section, and therefore in¬ 
dependent of the sign. As the magnitude is increased for 
the attractive potential, the velocity at which the dis¬ 
sipation increases steeply is only slightly reduced. But 
for the repulsive stirring potential a visible reduction is 
achieved, indicating the appearance of an additional dis¬ 
sipative mechanism. 

To understand the additional decay mechanism for the 
repulsive stirring potential, we show the time evolution 
of the phase field ^(r, t) of a single realization of the 
ensemble that is used in the simulation. This field is de¬ 
rived from the complex field '!/’('') via the phase-density 
representation ip = -^/n exp{i(p). In Fig. we show the 
time evolution of the phase field for a single ccy-plane 
of the three dimensional lattice, for both an attractive 
and a repulsive stirring potential. The magnitude of the 
stirrer is Vq = ±0.70 g.. In both cases, the velocity of 
the stirrer is above the steep onset of dissipation asso¬ 
ciated with the breakdown of superfluidity. For the re¬ 
pulsive stirrer we use v = 3.7 mm/s, for the attractive 
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FIG. 5. (Color online) The phase evolution of a single realization for a repulsive stirrer (upper row) and an attractive stirrer 
(lower row), at the times t = 26, 50, 71, and 195ms. The location of the stirrer is marked by a circle of radius a. For the 
repulsive stirrer the phase winding around the plaquettes is also indicated. A vortex (circled dot) and an anti-vortex (cross) 
correspond to a phase winding of -|-27r and —27r, respectively. For the repulsive stirrer, dissipation is caused by the creation of 
tightly bound vortex pairs, for the attractive stirrer by the creation of phonons. 


'O 






V > Vc 




vO 




FIG. 6. (Color online) Illustration of the dissipation mecha¬ 
nism of a repulsive stirrer, due to the creation of vortex-anti- 
vortex pairs. The stirrer is depicted by the blue disc, the 
arrow indicates the direction of the moving stirrer. A vortex 
(an anti-vortex) is depicted by the anticlockwise (clockwise) 
open circle arrow. 


stirrer we use v = 4.4mm/s. The time evolution of the 
phase during stirring displays strong phase gradients for 
repulsive stirring, while the attractive stirrer only cre¬ 
ates a smoothly varying phase field, associated with the 
creation of phonons. The phase evolution for the repul¬ 
sive stirrer is characterized by the creation of vortex- 
anti-vortex pairs. This can be confirmed by displaying 
the phase winding around each plaquette, in particular 
Ln H{x,y) = Sa:(l){x,y)+6y(j){x + l,y)+5a;(l)ix + l,y + l) + 
Sy(l){x,y + l), where the phase differences between sites is 
taken to be 5x/y(j){x,y) S (—7r,7r]. We display the calcu¬ 
lated phase winding in Fig. where closely bound pairs 
of phase winding -|- 27 r and — 27 r are observed, correspond¬ 
ing to vortex pairs. To display the decay with clarity, we 
also sketch it in Fig. We note that this mechanism 
of vortex-pair-induced dissipation has not been consid- 



FIG. 7. (Golor online) Simulated heating rates obtained 
for linear and circular stirring, with an attractive stirrer of 
strength Vo = —0.03 y. The onset of heating occurs at a crit¬ 
ical velocity Vc = 4.2 and 3.3mm/s for the linear and circular 
stirring, respectively. Here, Vc is determined using Eq. and 
the fitted curves are shown by the solid lines. 


ered or suggested in the literature. However, a recent 
experimental study displays similarities. Ref. |19j . 


V. CIRCULAR VERSUS LINEAR STIRRING 


In a trapped condensate, a circularly moving stirrer is 
a natural choice to probe the system at a fixed density. 
However, here we point out that the circular stirring mo¬ 
tion results in an intrinsic reduction of the onset of dis¬ 
sipation to lower velocities, compared to a linear stirring 
motion of the same velocity. In Fig. [^we show a compar¬ 
ison of the heating rate due to linear and circular stirring 
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FIG. 8. (Color online) (a) The temperature dependent critical velocities Vc (crosses) obtained from the simulation are compared 
to the critical velocities Vc,p (plus symbols) obtained from the analytical heating rate, Eq. [fusing the fitting function in Eq. 
1^ They are in good agreement but significantly below the sound velocity obtained from the simulation Vs (open circles) and 
from the condensate density at finite T, vb,t (filled squares). The gray, vertical thick line corresponds to the mean-field energy 
~ 0.56 ksTc- (b) The propagation of the outgoing density waves after the perturbation is turned off at time t = 0. (c) The 
sound velocity obtained in the simulation is influenced by the initial potential strength Vo. As is shown here for the various 
potential strengths To used in the simulation, and is compared to the sound velocity obtained from the local density excess 
WB.An (plus symbols) and from the applied potential ub.Vq (dashed line). 


motion. We calculate the heating rates of a homogeneous 
condensate at T = 0.07Tc. The circular motion of the 
stirrer has a radius i?o and a frequency Wm, the velocity 
of the stirrer is therefore given by u = RocOm- In this 
example, we choose Rq = 10 /tm and calculate the heat¬ 
ing rates as a function of stirring velocity v. The onset 
of dissipation for the circular stirring motion is reduced, 
compared to the linear motion. 

We now estimate the heating rate analytically, within 
Bogoliubov theory, in analogy to the previous results for 
a linear stirring motion. We consider a stirring potential 
of the form: 

= (23) 

with ys{t)) =_Ro(cos(wmt),sin(wmt)). As de¬ 

scribed in Appendix |B 3[ the equation of motion is now 
solved by 


This results in the heating rate 

jp A poo 

^ / dkkuJk{uk+Vk)'^No\Vk\'^J^{kRo) 

X S{uJk - (27) 

We approximate the sum over the Bessel function index 
v by an integral over a continuous variable, and obtain 

^ J dfcfcw/c(ufc-bUfc)^Afo|t4pT^(fci?o)- (28) 

To obtain the velocity scale of the onset of dissipation we 
note that for large indices (kRo) has its maximum 

at pT] 

kmRo- —+Co(—y^\ (29) 


= + (24) 

with 

2' ^ 

A'^,{t) = -jr{uk+Vk)^/^oVi, Y, i''e-^''^MkRo) 

iy— — (x> 

^ sin[iu:k - iyuJm)t/2] ^ 

{U!k - IXUJm) 


where J^{kRo) are the Bessel functions of the first kind 
of order ly. With this, and by integrating out the polar 
angle, the change of energy is 


9 A 

(Ho(t))-(Ho(0)) = — 


dk kujkiuk +Vk)‘^No\Vk 




sin^ [{uk - vuJm)t/2] 
{uJk - 


( 26 ) 


where Co « 0.809. Note that fci?o = is the same 

statement as vk = uik, i.e. we recover approximately the 
Landau criterion. Eq. can also be written as 

vk^ Ki u}k(l + (30) 


where Uc = v^/Rq is the centripetal acceleration. The 
maximum thus approaches the location of the 5-function 
in Eq. as Rq is increased. The width of that maxi¬ 
mum is tJB = 1/{2\/Cq) X {ookliyJm )^/^. The ratio of the 
width to the location of the maximum is 


ctb _ 1 / Oe \ 

kmRo 2 ^/(Tq V ) 


(31) 


This is the intrinsic broadening of the peak due to the 
accelerated nature of the circular motion. As we increase 















the radius Rq of the circular motion, the ratio of the 
centripetal acceleration and the radius adRo approaches 
zero. Thus, the function (kRo) becomes more and 

more narrow around its maximum, and the sharp onset 
of dissipation for a linearly moving stirrer is recovered. 

Using the expressions for Uk and (uk + Vk) in Eq. 28 
we obtain 


Rsc = 


27r^(T^ 


mAuj, 


.4 poo 

— / dkdexp{-k'^d)J^,/^^ikRo). (32) 

'm JQ 


We show this analytical prediction in Fig. ib) in com¬ 
parison to heating rates that were obtained numerically, 
which gives good agreement. We also show the analyt¬ 
ical prediction for th e latt ice system b y the dotted line 
in Fig. Hb), see Eq. |B48| in Appendix |B 3[ giving very 
good agreement. 

Written in dimensionless form, in analogy to Eq. 18 
we rewrite Eq. [3^ as 


i?n = VSttv 


pcrku 


d. - dkR^ja) 


(33) 


where k = ka, and v = Rd’^'eliyad- We show the 
comparison of the heating rates i?n for the linear and 
circular stirring in Fig. Both this figure and Fig. 
display the difference between the cases of linear and cir¬ 
cular stirring. For circular stirring, both the width and 
the location of the peak depend on the stirring radius 
i?o, described by Eq. Furthermore, the heating rate 
has a weak oscillatory behavior for large u, which is due 
to the oscillatory behavior of the Bessel functions, con¬ 
trolled by the radius i?o. As Rq is increased, the first peak 
builds up, and the oscillatory behavior shows a smaller 
and smaller period in u, visible in the inset of Fig. Htb). 
For very large Rq the onset of dissipation approaches that 
of linear stirring. 



FIG. 9. (Color online) The simulated heating rate (red 
squares connected with a dashed line) is compared to the 
analytical prediction (green, continuous lines) for linear stir¬ 
ring at finite temperatures: (a) T = 0.18 Tc, (b) T = 0.35 Tc, 
and (c) T = 0.57T,. The analytical results are obtained for a 
continuum system using Eq. |36| 


VI. INFLUENCE OF TEMPERATURE ON 

Up to now, we have considered temperatures that are 
small compared to the mean-field energy, and also small 
compared to the typical temperatures used in experi¬ 
ments, such as in Ref. [15]. 

In Fig. I^a) we show the temperature dependence of 
the critical velocity for a homogeneous condensate that 
is stirred with a linear motion. We use a lattice with 
128 X 32 X 16 sites and calculate the heating rates at dif¬ 
ferent temperatures T. From the heating rates we deter¬ 
mine the critical velocity using the fitting function in Eq. 

We use the linear stirrer with strength Vq = —0.03/r. 
We compare the temperature dependence of Vc to the Bo- 
goliubov velocity ub = ^/gnjm in Fig. ^A), where n is 
the total densit y. We c ompare Vc to theBogoliubov ve¬ 
locity vb,t = y/gnd^, where the numerically obtained 
ric is the condensate density at temperature T. At low 
temperatures and ub,t are close to each other, but 
start to deviate as the temperature is increased. As we 


demonstrate below, this is due to the thermal broadening 
of the phonon modes. 

In addition to the comparison to the estimated sound 
velocity we also determine the sound velocity Vs numer¬ 
ically. We use a Gaussian potential of width a along the 
cc-axis at the center of the homogeneous system, with 
a strength Vq = —0.03/i that we slowly ramp up over 
100 ms at the center and then abruptly turn it off. This 
results in a density wave traveling outward, as shown 
in Fig. [^b). The density wave splits into two outgo¬ 
ing waves propagating along the a:-axis. From the linear 
propagation we determine Vs- The example in Fig. ib) 
is for T = 0.07 Tc- We perform this simulation for a range 
of temperatures and show the temperature dependence 
in Fig. [^a). The results for Vs are in excellent agree¬ 
ment with vb,t- This demonstrates that the reduction of 
the critical velocity is not entirely due to the thermal re¬ 
duction of the sound velocity itself. We show below that 
the critical velocity can be recovered, if not only the re¬ 
duced sound velocity of the phonon modes but also their 
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thermal broadening is taken into account. 

We note that the velocities that are obtained in this 
way are influenced by the initial potential strength Vq. 
This is shown in Fig. |^c), at T = 0.07Tc, where we 
plot the velocity in units of wb as a function of Vq. As 
the strength Vb of the attractive potential is increased, 
the velocity Vs increases. For comparison we calculate 
the Bogoliubov velocity VB,An using the local density ex¬ 
cess An at the center owing to the ap plied potentia l Vq. 
In Fig. I^c) we show UB,vb/'*^B = Vl - |l^o|/(2/i) and 
VB,An/vB = -y/l -f An/(2n) by the dashed line and the 
plus symbols, respectively. Here, the factor 1/2 appears 
because the excited wave packet is divided into two equal 
density waves. For this reason, to obtain a reliable nu¬ 
merical estimate of the sound velocity, a sufficiently small 
value of the initial potential has to be used, as we did for 
the data depicted in Fig. |^a). 

To expand on the analytical approach given above, 
we include non-zero thermal broadening of the phonon 
modes due to Landau damping |21j . The ^-distribution 
In the Bogoliubov estimate of the heating rate, Eq. [T^ 
Is replaced by a Lorentzlan distribution with the width 
given by [2^ 


Ffc = X 


4 agkBTk 


h 


(34) 


Ffc is the width of a phonon mode with momentum k. 
The heating rate is then 


^ = ‘^^i^k{uk + Vk)^No\Vk 

at h 




(Wfc - vk)2 -H Tl 


■ (35) 


After substituting {uk + Vk)'^ and |I4;P, we rewrite Eq. 

155] as 


A,c = ^ / A^kk'^e 




(wfe - vkxY + r2 ’ 


(36) 


where = k^ + ky. We numerically evaluate Eq. 36 for 
a continuum system and compare this prediction with 
the simulated heating rate at different temperatures T in 
Fig. ^ The analytical results are In quantitative agree¬ 
ment with the simulations. To determine Vc we fit the 
analytical heating rate with Eq. We show the tem¬ 
perature dependence of the analytical in Fig. HJa). 
For temperatures T < /i/ZcB the predictions for Vc are in 
agreement with the simulated Vc- 


VII. UNIFORM VERSUS TRAPPED SYSTEM 


We now consider an inhomogeneous density distribu¬ 
tion, as it is realized in an atomic trap, and the influence 
of the inhomogeneity on the heating rate. As an exam¬ 
ple we show the numerically obtained heating rates for a 
trapped system in Fig. for parameters chosen in ac¬ 
cordance with the experiment |15j . corresponding to the 
experimental data point at —1/kpa « —3.5 in Fig. 11 



V [mm/s] 

FIG. 10. (Color online) The numerically obtained heating 
rates (circles), for a trapped condensate, are compared to the 
analytical predictions (dotted line). The analytical results use 
the effective condensate density obtained for a Thomas-Fermi 
profile. For comparison, we also show the numerical results 
for a homogeneous density of 0.42/rm“® and for 0.28 ^m“®, 
shown by the filled and empty squares, respectively. The den¬ 
sity of 0.42 /rm“® corresponds to the center density at stirrer 
location, the density of 0.28 /im“® corresponds to the effec¬ 
tive density at stirrer location. The inset shows the trapped 
and analytical results on a smaller scale, demonstrating good 
agreement at the onset of dissipation. 


The temperature is set to be T = 10 uK. For comparison, 
we also simulate a uniform system with an average den¬ 
sity of 0.42/im“^, which is the density uq at the stirrer 
location in the trap, see [23]. The onset of dissipation in 
the trapped system occurs at a smaller velocity as com¬ 
pared to the uniform case, see Fig. 

This leads us to consider an effective, reduced density 
instead of the maximal density of the trapped system. 
We determine the effective density Uq by averaging over 
a Thomas-Fermi distribution, (23], and obtain an effec¬ 
tive density of 2/3 times the central density. We now 
use this effective density in a simulation of a homoge¬ 
neous system. We show the results for the uniform case 
with the effective density ho at the stirrer location by 


the empty squares in Fig. 10 The comparison between 
the simulation for the trapped system and homogeneous 
system with the effective density shows good agreement. 

Furthermore, we use this effective density in the ana¬ 
lytical description of the heating process, in addition to 
extending the heating rate for circular stirring, Eq. |27| 
to finite temperatures by replacing the J-function with 
a Lorentzian distribution. Therefore, the finite tempera¬ 
ture heating rate is 


dE _ A ^ 


dfc kujkiuk + Vk)^ No\Vk\^ Jl{kR) 




(Wfc - ViOmY + Tfc ’ 

where the Landau damping F^ is given by Eq. 34 
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FIG. 11. (Color online) The experimentally measured (filled 
circles), the numerically simulated (crosses) and the analyt¬ 
ical prediction (plus symbols) for the critical velocity Vc in 
units of the Fermi velocity u_f in the BEG regime. The nu¬ 
merical results and the experimental results show quantita¬ 
tive agreement throughout the BEG regime, excluding the 
strongly correlated regime. The analytical prediction shows 
good agreement for the weak stirrer regime. In addition, we 
show the numerically obtained Vc for a weak stirrer of strength 
Vb = —h X 2.8 feHz, shown by the empty circle. This recovers 
the analytical prediction, which assumes a weak stirrer. 


Rewriting Eq. as 


R 


SC 



fc^exp(-fc^g^)jg(fcR)rfc 
{uJk - + rl 


We evaluate Eq. using the effective density h-o and 
the temperature T = 10 nK. We show this analytical 
prediction by the dotted line in Fig. [T^ This indeed 
captures the onset of dissipation well. 


VIII. COMPARISON TO EXPERIMENT 


Finally, we compare the simulation and the analyti¬ 
cal results to the experimental results [13] in the BEC 
regime. We use the same parameters as the experi¬ 
ment and determine the critical velocities Vc via the 
fitting method described before, for various scattering 
lengths tts = aoD, az>z> being the dimer-dimer scattering 
length. The dimer-dimer scattering length is calculated 
from the atom-atom scattering length a using the relation 
o-DD = 0.6a [531153]. The experimentally measured and 
the simulated critical velocities Vc in units of the Fermi 
velocity vp are plotted against the interaction strength 
in Fig. El The Fermi velocity of a non-interacting gas 
is given by vp = x rua is the 

mass of a ®Li atom, N is the total number of atoms, 
and uJr (w^) is the trap frequency in the radial (trans¬ 
verse) direction. The interaction strength is shown in 
terms of the dimensionless quantity —1/kpa, where a is 


the atom-atom scattering length and the Fermi wavevec- 
tor kp = niaVp/h. Throughout the BEC regime the 
experimentally measured and the numerically obtained 
Vc are in good agreement. We note that the simula¬ 
tions start to deviate from the experimental results in 
the strongly interacting limit —1/kpa > —1, where the 
simulation method becomes unreliable due to strong cor¬ 
relations. Furthermore, we display the analytical predic¬ 
tion based on the finite-temperature heating rate, Eq. 
|38[ for a circularly moving stirring potential, and for the 
effective density based on the Thomas-Fermi approxima¬ 
tion. From this heating rate we determine an estimate of 
the critical velocity, again via fitting. 

Again, we find remarkable agreement for interactions 
— 1 > —\/kpa > —2.8. The theoretical prediction starts 
to deviate both for the strongly-correlated regime, and 
for —1/kpa < —2.8. While the breakdown for the 
strongly-correlated regime is again due to the inapplica¬ 
bility of a weak-coupling method, the disagreement for 
weak interactions is due to entering the strong-stirrer 
regime, i.e. |Vb| > Mo is the mean-field energy at the 
stirrer location. The regime where the stirrer strength 
is comparable to the mean-field energy is shown by the 
gray, vertical thick line in Fig. 11 to the left of this 
line we are in the strong-stirrer limit. We demonstrate 
that this is the correct interpretation, by showing a weak 
stirrer given by Vq = —h x 2.8/cHz, in contrast to the 
experimentally used Vq = —h x 8.5fcHz, by the empty 
circle in Fig. E This simulation indeed reproduces the 
analytical prediction. 


IX. CONCLUSIONS 

In conclusion, we have demonstrated a quantitative 
understanding of the stability of superfluidity of a Bose- 
Einstein condensate stirred with an attractive potential. 
While we show that a similar, repulsive perturbation re¬ 
sults in a decay mechanism due to vortex-anti-vortex 
pairs, which had not been considered before, we And 
that the choice of an attractive perturbation results in 
phononic decay only, while the influence of vortices is 
suppressed. We then continue to identify the reason 
why, despite this phonon-only mechanism, the critical 
velocity that was obtained in the experiment m is still 
noticeably smaller than the phonon velocity itself. The 
contributing factors are the circular motion of the stir¬ 
rer, which results in an intrinsic broadening due to an 
accelerated motion, and the thermal broadening. With 
these factors taken properly into account, a quantitative, 
even analytical, agreement is achieved. The analytical 
prediction for the heating rate as a function of the stir¬ 
ring velocity is summarized in Fig. |12[ While the zero- 
temperature, linear stirring process gives the phonon ve¬ 
locity as the critical velocity, the onset of dissipation is 
noticeably reduced due to the thermal broadening and us¬ 
ing an accelerated stirrer. Beyond giving a quantitative 
understanding of the experiment reported in Ref. m, 
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FIG. 12. (Color online) Analytical heating rates Rn as a func¬ 
tion of stirring velocity v at temperature T = 0 and 10 nK. 
Panel (a) corresponds to linear stirring, whereas panel (b) 
shows the results for circular stirring with radius Ro = 10 /rm. 

this study thus provides a general blueprint to analyze 
the breakdown of superfluidity, disturbed with a weak, 
localized pertubation, in a typical experimental setting. 



FIG. 13. (Color online) The total energy change AE (blue 
circles) during stirring of a homogeneous system with a lin¬ 
ear stirrer of strength Vo = — 0.03 /r. The simulated energy 
is normalized by the total number of atoms. Panel (a) corre¬ 
sponds to a stirring velocity w < ub , whereas panel (b) shows 
the total energy change for v > vb- We fit the linear energy 
increase during stirring with a linear function. The fitting 
duration is marked by the two vertical, dashed lines, and the 
fitted curves are shown by the continuous lines. 
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Appendix A: Simulated heating rate 

In this section we illustrate how we obtain the heat¬ 
ing rate from the numerical simulation method described 
in Sec. [H] To perform the simulation we discretize the 
real space and represent the continuous Hamiltonian, Eq. 

by the discrete Bose-Hubbard Hamiltonian [57] on a 
three-dimensional square lattice: 

Ho = -j 

(b) » * 

where and = IV'iP are the complex-valued field and 
the density at site i, respectively; (ij) represents nearest- 
neighbor bonds. For a lattice discretization length I, the 
Bose-Hubbard parameters are related to the continuum 
parameters, cp. Ref. [inj, by J = h'^/{2ml'^) and U = 
gl~^. m is the atomic mass and g = jm is the 


interaction parameter; Og is the s-wave scattering length. 
The external potential Vi describes the harmonic trap, 
Vtrap.j = rn{ui'^r'^ +U!'^z'^)/2. ujr {uJz) is the trap frequency 
in the radial (transverse) direction, and r = {x^ V 
is the radial coordinate. 

As an example we choose the homogeneous system de¬ 
scribed in Sec. m with linear stirring. We initialize the 
system in a thermal state at temperature T. The initial 
states are generated from a grand canonical ensemble via 
a classical Metropolis algorithm. We then introduce the 
stirring potential described by Eq. [fusing classical equa¬ 
tions of motion. After following the stirring sequence de¬ 
scribed in Sec. m we calculate a thermal ensemble of the 
total energy E = (Hq) using the unperturbed Hamilto¬ 
nian in Eq. as a function of the stirring velocity v. 
We show two such cases, below and above the Bogoliubov 
velocity ub, for the linear stirring in Eig. |13[ During the 
stirring the energy increases linearly. Erom the linear en¬ 
ergy increase we determine the induced heating by the 
stirrer. 


Appendix B: Analytical heating rate 

We elaborate on calculations of the heating rate within 
the Bogoliubov approximation for both linear and cir¬ 
cular stirring. In Sec. |B 1| we describe the Bogoliubov 
Hamiltonian at temperature T = 0. In Sec. |B 2| we de- 
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rive the heating rate for linear stirring, Eqs. and 
In Sec. |B 3| we show the full calculations of the heating 
rate for circular stirring, Eqs. [27|and[^ 


Substituting V (r, t) using Eq. B6 and making use of the 
property LSk = f dr exp(—ikr), with L being the system 
length, we obtain 


1. Bogoliubov Hamiltonian 


Vk(t) 


27rVbcr^ ^ 

4. 


^iikxV:^+kyVy)t^-kl.^a^/2 


(B8) 


In this section we describe the Bogoliubov Hamilto¬ 
nian, the unperturbed Hamiltonian at temperature T = 
0, which is the starting point of our calculations. For a 
uniform gas of N interacting bosons in a box of volume 
V, the Hamiltonian in momentum space is given as 

no = Y^ efeoj^ak + ^ 4-(-p4-p“q“k, (BI) 

k k.q.p 

where Ok (aj;.) is the annihilation (creation) operator and 
ek = /2m. We diagonalize the Hamiltonian in Eq. 


UfcOk — ffcoLki where the Bogoliubov functions are given 
by ul = i^k + ek+ gric)/{2hujk) and vl = {-huk + ek + 
gUc) /(2IiWfe), respectively. The diagonalized Hamiltonian 
is therefore 


BI 


using the Bogoliubov transformation of the form 6k = 


Ho = hiOQ + ^ 6a;fc6k6k. 
k^O 


The dispersion of the quasiparticles is given by 


huJk = y e/c(efc + 2mu|), 
where vb = \/~gn//Jm is the sound velocity. 


(B2) 


(B3) 


2. Linear stirring 

Next, we introduce the stirring potential as a time- 
dependent perturbation. The total Hamiltonian of the 
system is then 

?i(t) =Ho + Hs(<), (B4) 

where the time-dependent Hamiltonian Us is given by 

Us{t) = j dvV{v,t)h(r). (B5) 

h(r) is the density operator and V (r, t) is the stirring 
potential. We model the stirrer as a Gaussian with a 
width a and a strength Vq. For a linear stirring motion 
the Gaussian stirrer moves with a velocity v = {vx,Vy). 
The stirring potential is therefore 

(B6) 

We Fourier transform the stirring potential using 

Vkit) = ^ J dre^’^Vir,t). (B7) 


where A is the system area and the Kronecker delta Sk^ 
is 1 for kz = 0 and zero otherwise. Rewriting Eq. |B8| as 

V^{t) = (B9) 

We now express the density operator h(r) in its Fourier 
representation, i.e. hk = X^k' ®k'®k'-i-k- Assuming a 
macroscopic occupation of the mode k = 0, we expand 

«k = aLk®o + ^ Syak'+k- (BIO) 

k'^O 


We linearize hk for small k using the Bogoliubov theory, 
where do and aj are replaced by y/No- Nq is the number 
of condensed atoms in the mode k = 0. Applying the 
Bogoliubov transformation, we obtain 

hk « \/f^(ufc-I-'(;fc)(&Lk + ^k). (BIl) 


Now, we rewrite the stirring potential and density in its 
Fourier representation in the stirring Hamiltonian as 

k k' 

= y drexp(-yk-tk')r). (BI2) 

k.k' 


\J Olilg, 


/ dr exp(—i(k -|- k')r), we get 


=5I^k(i)^ik = 5I'^'=(^-k + M. (BI3) 

k k 


Here, Sk = 27rVoi//Vo(uk -I- Vk)a^/A x 
The equation of motion for 6k (t) is written as 

indtb]^(t) = [6k(t),'Ho] + [6k(t),'Hs]. (BI4) 

We solve Eq. |B14| using the following ansatz: 

6k(t) = e-“''‘6k + Ak(t). (BI5) 

Using Eq. |B15| the terms on the right-hand side of Eq. 
IBI4l are 

[6k(t),'Ho] = [e"*‘^''‘6k-l-Ak(t),6a;o+ ^ hujk'bl,bi^>] 

k'^O 

= ^ huk' [6k, 6^6k-] = e-“'“‘6wfe6k, (BI6) 

k'^O 
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and 


[k{t),ns] = [e-^‘^'^% + Ak{t),J2Sk'{bl^,+k>)] 

k' 

J2Sk'[k,bl^^,] = e-^^'‘*Sk. (B17) 




Inserting Eqs. |B16 and B17 into Eq. B14[ which gives 
ihdtbi,{t) = e-*‘^’‘*huJkK + (B18) 

The time derivative of Eq. |B15| is 

ihdtbkit) = + iMtAk{t), 

which we equate with Eq. |B18| and get 


(B19) 


n 


= + Ufc) V^Vu. 


(B20) 


Here, Ek = ^-kV^o^/A x 5k^e~^ With the initial 

condition that the perturbation E(t < 0) = 0, for a finite 
stirring time t we obtain 


—2E , /— sin[(a;fe — vk)t/2l 

= —^{uk + Vk)y/^Vk - - -i 

h (ujk - vk) 


X e 


-i(uik-v'k)t /2 


(B21) 


We now calculate the energy of the system at time t 
using the dynamical evolution of the Bogoliubov operator 
created by Eqs. |B2| and [B5l The expectation value of the 
energy is 

(E(t))=^fiwfe(S|^(t)Sk(t)). (B22) 


Using the evolution of 6k (i) given in Eq. 


B15 


gives 




{E{t)) = Swfc((6j^6k) + Hk(t)(6j,)e' 

k 

+ H^(t)(6k)e-“'=* + |Hk(t)P). (B23) 


The first-order terms in Eq. B23 are zero, because (6k) = 
(^k) ~ d. Therefore, the energy change due to the second- 
order term in perturbation is 

{AE{t)) = y; HuJk\Ak{t)\'^. 


(B24) 


Substituting I Ak(t)r using Eq. B21 which gives 


(AU(t)) = -^Y ^k{uk + Ufc)^fVo|Uk 


sin^ [(wfc - vk)t/2] 

(Wfc - vk)2 


(B25) 


In the t 


oo limit we use 

sin^(at/2) tt 


lim 

t—^OO 




(B26) 


Using Eq. B26| the energy in Eq. B25 is proportional to 
t. Therefore, the heating rate is 

k 

Using the explicit form of |VkP a nd ( itfc -|- Vk)^ = 


B27 


we get 


hk^/{2mu!k), and J2k = 1 in Eq. 

4 p 

Rsc =/ d^fc exp(—A:^(T^)(5(a;fc — u/ca;), (B28) 

mA J 

where + ky, and R^c = {dE/dt)/{N qVq) is the 

rescaled heating rate. We have set the direction of the 
stirring velocity along the x-axis, i.e. v = ve^. In polar 
coordinates we integrate out the polar angle </>. This 
we do by replacing the d-function using limg_,,o = 

'kS{x). With this. 


Rkc. = 


27r(T^ 


exp(—fc^CT^) 

-1“ / dk — , 

mA Jq 


(B29) 


for vk > ojk, i.e. v > vb- The wavevector determines 
the maximum phonon momentum that can be created by 
a moving stirrer. Using the explicit form of ujk and in the 
ku ^ oo limit, we rewrite Eq. |B29|as 


,^'ECrik 


k^ exp(—fc^cr^) 


mA Jq — Ug — h?k‘^/{2m)‘^ 

Eq. |B30| can be solved analytically, giving 

Tr. 

Rue. = 


(B30) 


hA 


X Q(^-u(l/2,0,4rrPa‘^{vl - v‘^)/fp'^y (B31) 

where U{a,b,c) is confluent Hypergeometric function. 


3. Circular stirring 

In this section we derive the heating rate for a circular 
stirring motion with a radius Rq and a frequency 
The stirring potential has the following form: 

U(r,t) = Uoexp(-d- ^ ( 532 ) 

where (^Xs{t), ysit)) = Ro(cos{ujmt),sm{ujmt)) are the lo¬ 
cations of the Gaussian stirrer with width a and strength 
Vq. Eor a circular stirring motion, the Eourier transform 
of the potential is then 


Vi,{t) = A-t)+Kvs{t)) e-k^ A12 _ 


A 


(B33) 


X 
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In polar coordinates we obtain 

I4(t) = (B34) 

We now solve the equation of motion described by Eq. 
B14 using the stirring Hamiltonian defined by Eq. B5 


where we use the stirring potential for a circular stirring 
motion, Eq. |B32| The equation of motion is now solved 
by 




(B35) 


The solution of the equation of motion gives (similar to 
Eq. B20 for a linear stirring motion) 




_ _ —iujkt ikRQCOs{uJmt—(l)) 

~ H 
X {uk + 


(B36) 


Here, S'(, = [uk + Vk)\/^V\^ x exp cos{ojmt — </'))■ 
In Eq. |B36| we make use of the Jacobi Anger expansion: 


i’' . (B37) 


V — — CC) 


Ji,{x) is the Bessel function of the first kind of order v. 
After using Eq. B37 in Eq. B36 we get 


K{t) = -l{uk+ Vk)V^oVu Y J. ikRo) 




(B38) 


With the initial condition < 0) = 0 and for a finite 
stirring time t, this gives 

Kit) = 


sin[(a;fc - iyuJm)t/2] 
(Wfc - VUJm) 


. (B39) 


The energy of the s ystem due to the perturbation at 
time t is (see Eq. B22|) 


{AE{t)) = Y^>^\Kit)\^- 


(B40) 


Replacing |Aj^(t)p using Eq. B39 which gives 


{AE{t)) =-Y^kiuk+vKNo\Vi,\'^Y 

k I/',!/' 

sin[(a;fc - vuJm)t/2] sin[(a;fc - u'uim)t/2] 

(Wfc - lyUJm) (Wfc - v'uJm) 


(B4I) 


We use J2k = 1 in Eq. B41 and that reduces k = 
{kx,ky). In polar coordinates, we then integrate out the 


polar angle cj) using d(j)exp[—i{iy — = 2TrS,ji^^. 

With this, we arrive at 

{AE{t)) = — J dk kojkjuk +-Cfc)^A^o|Vkpy^./3(fcRo) 

sin^liujk — i^ujjn)t/2] 

X -^- J A (B42) 

[UJk - VOJm)^ 

In the t —> oo limit the steady value of the energy in Eq. 


B42 is proportional to t. Therefore, the heating rate is 
dk kujkiuk + Vk)'^No\Vk\‘^ J^{kRo) 


dE A ^ 


X 8{u)k - VUJm)- (B43) 

We solve the J-function in Eq. |B43 for variable k using 


the property S{f{x)) = 5{x — xq)/\ f'{xo)\. With this, we 
obtain 

^ = ^^koUJkoiuko+Vkoft^oKo^JUKRo) 


A/(2mwB)^ + (Sfco)2 


(B44) 


2(r?WB)^ + (fifco)^ 

with fco = (•\/4(muB)^ + {2mhvuJm)'^ — ‘2{mvB)‘^)^^‘^/fi- 
We now approximate the sum over the Bessel function 
index v by an integral over a continuous variable in Eq. 
|B43[ and then solve the J-function. This gives 

dE A 


dt tujj. 


dkkuJkiuk + Vk) No\V-i^ 


m JO 


X K/u:Jt^Ro)- 


(B45) 

We now express the heating rate for a lattice system. 
We therefore rewrite the potential in Eq. |B33| explicitly 
with kx, ky as 

Vk{t) = Ek 

X exp[iY^fc^~+~^i?o cos(wm^ — tan“^(A:j^/fca;))] . (B46) 
With this, the energy change in Eq. |B40| at time t is 

exp(-fc^tT^) Y JvikRo) 


{AE{t)) STr^cr"^ 




X J„>{kRo)cos {v - v')(J^+ - tan ^{ky/kx 


IX,I/' 

.-I/' 


sin[(a;fc - vuJm)tl2\ sin[(a;fc - v'uJm)t/2\ 


(B47) 


{uJk - VUJm) (Wfe - v'uJm) 

we used the explicit f orms of (uk + Vk)"^ and 


B47 


In Eq. 

iRkP- The time derivative of Eq. B47 gives the heating 
rate, which can be written as 

47r^ cr"^ 




mA'^ 


Y ^k^k"^ exp(-fc2cr^) Y J^if^Ro) 


sm\ (UJk — VUJm)t] , , , 

X -4^ -^ + Rsciv ^ v'). (B48) 


(Wfe - VUJm) 
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At long times, the heating rate is dominated by the di¬ 
agonal contributions, i.e. the terms with v = v'. We 


evaluate this part of the sum in Eq. |B48| for the first 
Brillouin zone, and show this analytical prediction by 
the dotted line in Fig. ib). 
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